Intro

In this script, I fit simplified models to density data so that I can predict densities on the condition data and pred grid

Load libraries

rm(list = ls())

library(tidyverse); theme_set(theme_light(base_size = 12))
library(readxl)
library(tidylog)
library(RCurl)
library(viridis)
library(RColorBrewer)
library(patchwork)
library(janitor)
library(icesDatras)
library(mapdata)
library(patchwork)
library(rgdal)
library(raster)
library(sf)
library(rgeos)
library(chron)
library(lattice)
library(ncdf4)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(geosphere)
#remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB)
# To load entire cache in interactive r session, do:
# qwraps2::lazyload_cache_dir(path = "R/clean_data/cod_fle_density_models_as_covars_cache/html")

For maps

world <- ne_countries(scale = "medium", returnclass = "sf")

# Specify map ranges
ymin = 54; ymax = 58; xmin = 12; xmax = 22

map_data <- rnaturalearth::ne_countries(
  scale = "medium",
  returnclass = "sf", continent = "europe")

# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
  st_crop(map_data,
          c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))

# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)

#ggplot(swe_coast_proj) + geom_sf()

Load data

d <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cpue_q_1_4.csv")

# Filter to match the data I want to predict on
d <- d %>% filter(year > 1992 & year < 2020 & quarter == 4)

# Calculate standardized variables
d <- d %>% 
  mutate(depth_sc = (depth - mean(depth))/sd(depth)) %>%
  mutate(year = as.integer(year)) %>% 
  drop_na(depth) %>% 
  rename("density_cod" = "density") # to fit better with how flounder is named

# Explore data a bit
ggplot(d) + 
  geom_point(aes(year, density_cod, color = factor(sub_div))) + 
  facet_wrap(~sub_div)


# Explore data a bit
ggplot(d) + 
  geom_point(aes(year, density_fle, color = factor(sub_div))) + 
  facet_wrap(~sub_div)

Load prediction grid

pred_grid1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(1_2).csv")
#> Rows: 129346 Columns: 29
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (1): ices_rect
#> dbl (28): X, Y, depth, year, oxy, temp, lon, lat, sub_div, density_saduria, ...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pred_grid <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(2_2).csv")
#> Rows: 120107 Columns: 29
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (1): ices_rect
#> dbl (28): X, Y, depth, year, oxy, temp, lon, lat, sub_div, density_saduria, ...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

pred_grid <- bind_rows(pred_grid1, pred_grid)

# Standardize data with respect to data
pred_grid <- pred_grid %>%
  mutate(year = as.integer(year)) %>% 
  filter(year %in% c(unique(d$year))) %>% 
  mutate(depth_sc = (depth - mean(d$depth))/sd(d$depth),
         temp_sc = (temp - mean(d$temp))/sd(d$temp),
         oxy_sc = (oxy - mean(d$oxy))/sd(d$oxy)) %>% # Need to scale these to the mean and sd in the data!
  drop_na(oxy, depth, temp)
#> mutate: converted 'year' from double to integer (0 new NA)
#> filter: no rows removed
#> mutate: new variable 'depth_sc' (double) with 6,238 unique values and 0% NA
#>         new variable 'temp_sc' (double) with 249,453 unique values and 0% NA
#>         new variable 'oxy_sc' (double) with 249,453 unique values and 0% NA
#> drop_na: no rows removed

Make mesh

# Make mesh:
spde <- make_mesh(d, xy_cols = c("X", "Y"),
                  n_knots = 200, 
                  type = "kmeans", seed = 42)

plot(spde)

Fit models

# Depth spline + oxy spline
# Cod model
mcod <- sdmTMB(density_cod ~ 0 + as.factor(year) + s(depth_sc),
               data = d,
               mesh = spde, family = tweedie(link = "log"),
               spatiotemporal = "AR1",
               spatial = "on",
               time = "year",
               reml = FALSE,
               control = sdmTMBcontrol(newton_steps = 1))

tidy(mcod, conf.int = TRUE)
#>                   term estimate std.error conf.low conf.high
#> 1  as.factor(year)1993 4.837552 0.4688664 3.918590  5.756513
#> 2  as.factor(year)1994 5.197526 0.4517869 4.312040  6.083012
#> 3  as.factor(year)1995 5.807254 0.4443790 4.936288  6.678221
#> 4  as.factor(year)1996 5.083857 0.4449677 4.211736  5.955977
#> 5  as.factor(year)1997 4.316940 0.4287832 3.476540  5.157339
#> 6  as.factor(year)1998 4.502834 0.4325650 3.655022  5.350645
#> 7  as.factor(year)1999 4.347810 0.4224004 3.519921  5.175700
#> 8  as.factor(year)2000 4.694878 0.4021157 3.906746  5.483011
#> 9  as.factor(year)2001 5.099759 0.3908475 4.333712  5.865806
#> 10 as.factor(year)2002 5.712133 0.3823567 4.962728  6.461539
#> 11 as.factor(year)2003 4.502260 0.4054629 3.707568  5.296953
#> 12 as.factor(year)2004 4.959874 0.4083146 4.159592  5.760156
#> 13 as.factor(year)2005 5.346098 0.3743546 4.612376  6.079819
#> 14 as.factor(year)2006 5.041373 0.3600782 4.335632  5.747113
#> 15 as.factor(year)2007 5.280502 0.3580791 4.578680  5.982324
#> 16 as.factor(year)2008 5.495251 0.3536101 4.802188  6.188314
#> 17 as.factor(year)2009 5.438676 0.3608035 4.731514  6.145838
#> 18 as.factor(year)2010 5.418556 0.3576511 4.717572  6.119539
#> 19 as.factor(year)2011 4.679926 0.3587164 3.976855  5.382997
#> 20 as.factor(year)2012 4.422424 0.3662989 3.704492  5.140357
#> 21 as.factor(year)2013 4.735018 0.3622512 4.025019  5.445017
#> 22 as.factor(year)2014 4.582177 0.3638320 3.869079  5.295275
#> 23 as.factor(year)2015 4.681539 0.3634106 3.969268  5.393811
#> 24 as.factor(year)2016 4.124882 0.3637063 3.412030  4.837733
#> 25 as.factor(year)2017 4.560351 0.3605213 3.853742  5.266960
#> 26 as.factor(year)2018 3.664196 0.3665505 2.945770  4.382621
#> 27 as.factor(year)2019 3.762937 0.4148706 2.949806  4.576069
sanity(mcod)
#> ✓ Non-linear minimizer suggests successful convergence
#> ✓ Hessian matrix is positive definite
#> ✓ No extreme or very small eigen values detected
#> x `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> x `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> x `ln_tau_O` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> x `ln_tau_E` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> x `ln_kappa` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> x `ln_smooth_sigma` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✓ No fixed-effect standard errors are NA
#> ✓ No fixed-effect standard errors look unreasonably large
#> ✓ No sigma parameters are < 0.01
#> ✓ No sigma parameters are > 100
#> ✓ Range parameter doesn't look unreasonably large
mcod_1 <- run_extra_optimization(mcod, nlminb_loops = 1, newton_loops = 1)

sanity(mcod_1)
#> ✓ Non-linear minimizer suggests successful convergence
#> ✓ Hessian matrix is positive definite
#> ✓ No extreme or very small eigen values detected
#> ✓ No gradients with respect to fixed effects are >= 0.001
#> ✓ No fixed-effect standard errors are NA
#> ✓ No fixed-effect standard errors look unreasonably large
#> ✓ No sigma parameters are < 0.01
#> ✓ No sigma parameters are > 100
#> ✓ Range parameter doesn't look unreasonably large
mcmc_res <- residuals(mcod, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200)
#> Constructing atomic tweedie_logW
#> Constructing atomic invpd
#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.006393 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 63.93 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 1041.76 seconds (Warm-up)
#> Chain 1:                6.70257 seconds (Sampling)
#> Chain 1:                1048.47 seconds (Total)
#> Chain 1:

qqnorm(mcmc_res);qqline(mcmc_res)

# Save model (so that I can add density predictions to the pred grid for condition)
saveRDS(mcod_1, "output/mcod.rds")
# Fit flounder model
mfle <- sdmTMB(density_fle ~ 0 + as.factor(year) + s(depth_sc),
               data = d,
               mesh = spde,
               family = tweedie(link = "log"),
               spatiotemporal = "AR1",
               spatial = "on",
               time = "year",
               reml = FALSE,
               control = sdmTMBcontrol(newton_steps = 1))

tidy(mfle, conf.int = TRUE)
#>                   term estimate std.error conf.low conf.high
#> 1  as.factor(year)1993 3.556486 0.6490165 2.284437  4.828535
#> 2  as.factor(year)1994 3.083573 0.6494050 1.810763  4.356384
#> 3  as.factor(year)1995 4.232611 0.6346263 2.988767  5.476456
#> 4  as.factor(year)1996 3.618856 0.6370869 2.370188  4.867523
#> 5  as.factor(year)1997 3.397782 0.6169767 2.188530  4.607034
#> 6  as.factor(year)1998 2.783891 0.6189851 1.570703  3.997080
#> 7  as.factor(year)1999 2.350747 0.6129248 1.149437  3.552058
#> 8  as.factor(year)2000 3.154737 0.5957873 1.987016  4.322459
#> 9  as.factor(year)2001 3.575064 0.5883076 2.422002  4.728125
#> 10 as.factor(year)2002 4.537331 0.5806431 3.399291  5.675371
#> 11 as.factor(year)2003 3.375738 0.5913574 2.216699  4.534777
#> 12 as.factor(year)2004 3.844600 0.5955993 2.677247  5.011953
#> 13 as.factor(year)2005 3.704495 0.5725463 2.582324  4.826665
#> 14 as.factor(year)2006 3.668069 0.5594159 2.571634  4.764504
#> 15 as.factor(year)2007 3.942359 0.5580117 2.848676  5.036041
#> 16 as.factor(year)2008 4.099267 0.5556642 3.010185  5.188349
#> 17 as.factor(year)2009 4.186882 0.5603307 3.088654  5.285110
#> 18 as.factor(year)2010 4.225439 0.5552249 3.137218  5.313659
#> 19 as.factor(year)2011 3.853699 0.5550766 2.765769  4.941629
#> 20 as.factor(year)2012 3.604305 0.5573405 2.511938  4.696673
#> 21 as.factor(year)2013 3.837236 0.5587290 2.742147  4.932324
#> 22 as.factor(year)2014 3.606003 0.5583835 2.511592  4.700415
#> 23 as.factor(year)2015 3.746447 0.5577288 2.653319  4.839575
#> 24 as.factor(year)2016 3.768052 0.5565157 2.677301  4.858803
#> 25 as.factor(year)2017 3.901198 0.5575261 2.808466  4.993929
#> 26 as.factor(year)2018 3.636791 0.5605948 2.538045  4.735536
#> 27 as.factor(year)2019 3.920627 0.5897827 2.764674  5.076580
sanity(mfle)
#> ✓ Non-linear minimizer suggests successful convergence
#> ✓ Hessian matrix is positive definite
#> ✓ No extreme or very small eigen values detected
#> x `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> x `b_j` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> x `ln_phi` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✓ No fixed-effect standard errors are NA
#> ✓ No fixed-effect standard errors look unreasonably large
#> ✓ No sigma parameters are < 0.01
#> ✓ No sigma parameters are > 100
#> ✓ Range parameter doesn't look unreasonably large
mfle_1 <- run_extra_optimization(mfle, nlminb_loops = 1, newton_loops = 1)
sanity(mfle_1)
#> ✓ Non-linear minimizer suggests successful convergence
#> ✓ Hessian matrix is positive definite
#> ✓ No extreme or very small eigen values detected
#> ✓ No gradients with respect to fixed effects are >= 0.001
#> ✓ No fixed-effect standard errors are NA
#> ✓ No fixed-effect standard errors look unreasonably large
#> ✓ No sigma parameters are < 0.01
#> ✓ No sigma parameters are > 100
#> ✓ Range parameter doesn't look unreasonably large
mcmc_res_fle <- residuals(mfle_1, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200)
#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.006621 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 66.21 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 1030.14 seconds (Warm-up)
#> Chain 1:                5.2098 seconds (Sampling)
#> Chain 1:                1035.35 seconds (Total)
#> Chain 1:
#> Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
#> http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
#> Warning: Examine the pairs() plot to diagnose sampling problems

qqnorm(mcmc_res_fle);qqline(mcmc_res_fle)

# Save model (so that I can add density predictions to the pred grid for condition)
saveRDS(mfle_1, "output/mfle.rds")

Now that the models have been fit, I can create the pred_gid (make_pred_grid_utm.R)

[Extra for bentfish - calculate biomass index per subdivision]

### Cod
# SD 24
preds_cod24_sim <- predict(mcod_1, newdata = filter(pred_grid, sub_div == 24), nsim = 100)
#> filter: removed 211,032 rows (85%), 38,421 rows remaining

# SD 25
preds_cod25_sim <- predict(mcod_1, newdata = filter(pred_grid, sub_div == 25), nsim = 100)
#> filter: removed 178,983 rows (72%), 70,470 rows remaining

# SD 26
preds_cod26_sim <- predict(mcod_1, newdata = filter(pred_grid, sub_div == 26), nsim = 100)
#> filter: removed 184,410 rows (74%), 65,043 rows remaining

# SD 27
preds_cod27_sim <- predict(mcod_1, newdata = filter(pred_grid, sub_div == 27), nsim = 100)
#> filter: removed 224,505 rows (90%), 24,948 rows remaining

# SD 28
preds_cod28_sim <- predict(mcod_1, newdata = filter(pred_grid, sub_div == 28), nsim = 100)
#> filter: removed 198,882 rows (80%), 50,571 rows remaining

# Now calculate the CPUE index (average)
index24_cod_sim <- get_index_sims(preds_cod24_sim, area = rep(2*2, nrow(preds_cod24_sim)))
index25_cod_sim <- get_index_sims(preds_cod25_sim, area = rep(2*2, nrow(preds_cod25_sim)))
index26_cod_sim <- get_index_sims(preds_cod26_sim, area = rep(2*2, nrow(preds_cod26_sim)))
index27_cod_sim <- get_index_sims(preds_cod27_sim, area = rep(2*2, nrow(preds_cod27_sim)))
index28_cod_sim <- get_index_sims(preds_cod28_sim, area = rep(2*2, nrow(preds_cod28_sim)))

# Add some columns
index24_cod_sim <- index24_cod_sim %>% mutate(sub_div = "24", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index25_cod_sim <- index25_cod_sim %>% mutate(sub_div = "25", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index26_cod_sim <- index26_cod_sim %>% mutate(sub_div = "26", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index27_cod_sim <- index27_cod_sim %>% mutate(sub_div = "27", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index28_cod_sim <- index28_cod_sim %>% mutate(sub_div = "28", species = "cod")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA


### Flounder
# SD 24
preds_fle24_sim <- predict(mfle_1, newdata = filter(pred_grid, sub_div == 24), nsim = 100)
#> filter: removed 211,032 rows (85%), 38,421 rows remaining

# SD 25
preds_fle25_sim <- predict(mfle_1, newdata = filter(pred_grid, sub_div == 25), nsim = 100)
#> filter: removed 178,983 rows (72%), 70,470 rows remaining

# SD 26
preds_fle26_sim <- predict(mfle_1, newdata = filter(pred_grid, sub_div == 26), nsim = 100)
#> filter: removed 184,410 rows (74%), 65,043 rows remaining

# SD 27
preds_fle27_sim <- predict(mfle_1, newdata = filter(pred_grid, sub_div == 27), nsim = 100)
#> filter: removed 224,505 rows (90%), 24,948 rows remaining

# SD 28
preds_fle28_sim <- predict(mfle_1, newdata = filter(pred_grid, sub_div == 28), nsim = 100)
#> filter: removed 198,882 rows (80%), 50,571 rows remaining

# Now calculate the CPUE index (average)
index24_fle_sim <- get_index_sims(preds_fle24_sim, area = rep(2*2, nrow(preds_fle24_sim)))
index25_fle_sim <- get_index_sims(preds_fle25_sim, area = rep(2*2, nrow(preds_fle25_sim)))
index26_fle_sim <- get_index_sims(preds_fle26_sim, area = rep(2*2, nrow(preds_fle26_sim)))
index27_fle_sim <- get_index_sims(preds_fle27_sim, area = rep(2*2, nrow(preds_fle27_sim)))
index28_fle_sim <- get_index_sims(preds_fle28_sim, area = rep(2*2, nrow(preds_fle28_sim)))

# Add some columns
index24_fle_sim <- index24_fle_sim %>% mutate(sub_div = "24", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index25_fle_sim <- index25_fle_sim %>% mutate(sub_div = "25", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index26_fle_sim <- index26_fle_sim %>% mutate(sub_div = "26", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index27_fle_sim <- index27_fle_sim %>% mutate(sub_div = "27", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA
index28_fle_sim <- index28_fle_sim %>% mutate(sub_div = "28", species = "fle")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
#>         new variable 'species' (character) with one unique value and 0% NA

div_index_sim <- bind_rows(index24_cod_sim, index25_cod_sim, index26_cod_sim, index27_cod_sim, index28_cod_sim,
                           index24_fle_sim, index25_fle_sim, index26_fle_sim, index27_fle_sim, index28_fle_sim) %>% 
  mutate(est_t = est/1000, lwr_t = lwr/1000, upr_t = upr/1000) # convert to tonnes 
#> mutate: new variable 'est_t' (double) with 270 unique values and 0% NA
#>         new variable 'lwr_t' (double) with 270 unique values and 0% NA
#>         new variable 'upr_t' (double) with 270 unique values and 0% NA

Plot the index and save as csv

ggplot(div_index_sim, aes(year, est_t, color = sub_div)) +
  geom_line() +
  facet_wrap(~species, scales = "free") +
  geom_ribbon(aes(year, ymin = lwr_t, ymax = upr_t, fill = sub_div), alpha = 0.2) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
  guides(color = "none") + 
  labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "bottom",
        legend.key.height = unit(0.2, "cm"),
        legend.key.width = unit(0.2, "cm"),
        legend.background = element_blank(),
        plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
  NULL


ggplot(div_index_sim, aes(year, est_t, color = sub_div)) +
  geom_line() +
  facet_wrap(~species, scales = "free") +
  #geom_ribbon(aes(year, ymin = lwr_t, ymax = upr_t, fill = sub_div), alpha = 0.2) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
  guides(color = "none") + 
  labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "bottom",
        legend.key.height = unit(0.2, "cm"),
        legend.key.width = unit(0.2, "cm"),
        legend.background = element_blank(),
        plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
  NULL


ggplot(div_index_sim, aes(year, est_t, color = species)) +
  geom_line() +
  facet_wrap(~sub_div, scales = "free") +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2", name = "ICES subdivision") +
  labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
  theme(axis.text.x = element_text(angle = 90),
        legend.position = "bottom",
        legend.key.height = unit(0.2, "cm"),
        legend.key.width = unit(0.2, "cm"),
        legend.background = element_blank(),
        plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
  NULL


write.csv(div_index_sim, "output/cod_fle_index.csv")
knitr::knit_exit()